Hetero Analysis

n_distinct(result$MA.id) # number of MAs

## [1] 23492

nrow(result) # 9 model result per each MA

## [1] 211428

table((result %>% distinct(MA.id, .keep_all = TRUE))$event)

## 
##     common infrequent       rare  very rare 
##      22074       1154        101        163

proportions(table((result %>% distinct(MA.id, .keep_all = TRUE))$event)) * 100

## 
##     common infrequent       rare  very rare 
## 93.9639026  4.9123106  0.4299336  0.6938532

Estimation

Analysis of heterogeneity estimate

There are three types of model that have \(\tau\)

1. \(\tau\) from ivre, htevaguelogit and htelogit_tau50 (log ratios)

2. \(\tau_a\) and \(\tau_c\) from ablogit and ablogit_tau50 (log odds)

3. \(\tau_a\) and \(\tau_c\) from htebeta (prob)

Correlation of tau

cor(
  result_bayesian[result_bayesian$modelname == "htelogit_tau50",]$tau,
  result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$taua
  )

## [1] 0.2477361

cor(
  result_bayesian[result_bayesian$modelname == "htelogit_tau50",]$tau,
  result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$tauc
  )

## [1] 0.3046824

cor(
  result_bayesian[result_bayesian$modelname == "htebeta",]$taua,
  result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$taua
  )

## [1] -0.08585503

cor(
  result_bayesian[result_bayesian$modelname == "htebeta",]$tauc,
  result_bayesian[result_bayesian$modelname == "ablogit_tau50",]$tauc
  )

## [1] -0.1115232

Analysis Significance Inference by model

modelname	Sig_sum	percentage_total
ivfe	10403	44.283160
ivre	9132	38.872808
ctevaguelogit	11553	49.178444
ctebeta	11402	48.535672
htevaguelogit	7308	31.108462
htelogit_tau50	6047	25.740678
htebeta	4363	18.572280
ablogit	8003	34.066916
ablogit_tau50	597	2.541291

table(result_sigsum$Sig_sum)

## 
##     0     1     2     3     4     5     6     7     8     9 
## 10947  1075   706  1407  1587   922   915  2312  3121   500

round(proportions(table(result_sigsum$Sig_sum))*100,3)

## 
##      0      1      2      3      4      5      6      7      8      9 
## 46.599  4.576  3.005  5.989  6.755  3.925  3.895  9.842 13.285  2.128

table(result_sigsum$event, result_sigsum$Sig_sum)

##             
##                 0    1    2    3    4    5    6    7    8    9
##   common     9948  977  632 1340 1524  902  883 2261 3109  498
##   infrequent  799   74   58   59   61   18   29   42   12    2
##   rare         65    8    7    5    2    2    3    9    0    0
##   very rare   135   16    9    3    0    0    0    0    0    0

round(proportions(table(result_sigsum$event, result_sigsum$Sig_sum), margin = 1)*100,2)

##             
##                  0     1     2     3     4     5     6     7     8     9
##   common     45.07  4.43  2.86  6.07  6.90  4.09  4.00 10.24 14.08  2.26
##   infrequent 69.24  6.41  5.03  5.11  5.29  1.56  2.51  3.64  1.04  0.17
##   rare       64.36  7.92  6.93  4.95  1.98  1.98  2.97  8.91  0.00  0.00
##   very rare  82.82  9.82  5.52  1.84  0.00  0.00  0.00  0.00  0.00  0.00

waic

For best model comparison among seven Bayesian models (exclude ivfe and ivre)

htelogit

result_htevaguelogit_tau = result %>% filter(modelname %in% c("ctevaguelogit", "htevaguelogit")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>%  filter(modelname == "htevaguelogit") %>% ungroup()

mean(result_htevaguelogit_tau$diff_waic < 0) # htevaguelogit has better waic than ctevaguelogit

## [1] 0.9258471

result_htelogit_tau50_tau = result %>% filter(modelname %in% c("ctevaguelogit", "htelogit_tau50")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>%  filter(modelname == "htelogit_tau50") %>% ungroup()

mean(result_htelogit_tau50_tau$diff_waic < 0) # htelogit_tau50 has better waic than ctevaguelogit

## [1] 0.919164

result_htelogit = result %>% filter(modelname %in% c("htevaguelogit", "htelogit_tau50")) %>% group_by(MA.id) %>% mutate(diff_waic = diff(waic)) %>%  filter(modelname == "htelogit_tau50") %>% ungroup()

mean(result_htelogit$diff_waic < 0) # htelogit_tau50 has better waic than htevaguelogit

## [1] 0.5440576

ablogit

result_ablogit_tau = result_narm %>% filter(modelname %in% c("ctevaguelogit", "ablogit")) %>% group_by(MA.id) %>% 
  mutate(diff_waic = diff(waic)) %>%  filter(modelname == "ablogit") %>% ungroup()

mean(result_ablogit_tau$diff_waic < 0) # ablogit has better waic than ctevaguelogit

## [1] 0.5141309

result_ablogit_tau50_tau = result_narm %>% 
  filter(modelname %in% c("ctevaguelogit", "ablogit_tau50")) %>% group_by(MA.id) %>% 
  mutate(diff_waic = diff(waic)) %>%  filter(modelname == "ablogit_tau50") %>% ungroup()

mean(result_ablogit_tau50_tau$diff_waic < 0) # ablogit_tau50 has better waic than ctevaguelogit

## [1] 0.6156923

result_ablogit_tau = result_narm %>% filter(modelname %in% c("ctebeta", "ablogit")) %>% group_by(MA.id) %>% 
  mutate(diff_waic = diff(waic)) %>%  filter(modelname == "ablogit") %>% ungroup()

mean(result_ablogit_tau$diff_waic < 0) # ablogit has better waic than ctebeta

## [1] 0.7459887

result_ablogit_tau50_tau = result_narm %>% 
  filter(modelname %in% c("ctebeta", "ablogit_tau50")) %>% group_by(MA.id) %>% 
  mutate(diff_waic = diff(waic)) %>%  filter(modelname == "ablogit_tau50") %>% ungroup()

mean(result_ablogit_tau50_tau$diff_waic < 0) # ablogit_tau50 has better waic than ctebeta

## [1] 0.7681851

result_ablogit = result_narm %>% 
  filter(modelname %in% c("ablogit", "ablogit_tau50")) %>% group_by(MA.id) %>% 
  mutate(diff_waic = diff(waic)) %>%  filter(modelname == "ablogit_tau50") %>% ungroup()

mean(result_ablogit$diff_waic < 0) # ablogit_tau50 has better waic than ablogit

## [1] 0.6549662

##                ctevaguelogit ctebeta htevaguelogit ablogit htebeta
## ctevaguelogit          0.000   0.741         0.074   0.486   0.280
## ctebeta                0.259   0.000         0.062   0.254   0.040
## htevaguelogit          0.926   0.938         0.000   0.920   0.754
## ablogit                0.514   0.746         0.080   0.000   0.265
## htebeta                0.720   0.960         0.246   0.735   0.000
## htelogit_tau50         0.919   0.945         0.541   0.925   0.782
## ablogit_tau50          0.616   0.768         0.199   0.655   0.443
##                htelogit_tau50 ablogit_tau50
## ctevaguelogit           0.081         0.384
## ctebeta                 0.055         0.232
## htevaguelogit           0.459         0.801
## ablogit                 0.075         0.345
## htebeta                 0.218         0.557
## htelogit_tau50          0.000         0.828
## ablogit_tau50           0.172         0.000

##                ctevaguelogit ctebeta htevaguelogit ablogit htebeta
## ctevaguelogit          0.000   0.000         0.000   0.000   0.000
## ctebeta                0.259   0.000         0.000   0.000   0.000
## htevaguelogit          0.926   0.938         0.000   0.000   0.000
## ablogit                0.514   0.746         0.080   0.000   0.000
## htebeta                0.720   0.960         0.246   0.735   0.000
## htelogit_tau50         0.919   0.945         0.541   0.925   0.782
## ablogit_tau50          0.616   0.768         0.199   0.655   0.443
##                htelogit_tau50 ablogit_tau50
## ctevaguelogit           0.000             0
## ctebeta                 0.000             0
## htevaguelogit           0.000             0
## ablogit                 0.000             0
## htebeta                 0.000             0
## htelogit_tau50          0.000             0
## ablogit_tau50           0.172             0